Lab8_regressions

Регрессии для непрерывных зависимых переменных

Регрессия

Регрессия - метод предсказания (объяснения) одной переменной через другие.

\[y_i = \hat\beta_0 + \hat\beta_1 \times x_i + \epsilon_i,\] * независимые переменные (предикторы, effect) * зависимая переменная (outcome)

В отличие от корреляции, проверяющей связь между двумя переменными, но не направление зависимости, регрессия явно эксплицирует это направление.

library(tidyverse)
library(lme4)
library(lmerTest) # for the model summary
library(skimr) # just for summaries

set.seed(42)
data <- tibble(x = rnorm(150)+1) %>% 
  mutate(y = 5*x+10+rnorm(100, sd = 2))
mean_y <- mean(data$y)

data %>%
  ggplot(aes(x, y)) +
  geom_point() +
  geom_hline(yintercept = 0, linetype = 2) +
  geom_vline(xintercept = 0, linetype = 2) +
  geom_hline(yintercept = 15.314, color = "blue") +
  annotate(geom = "point", x = 0, y = 15.314, size = 4, color = "red") +
  scale_x_continuous(breaks = -2:4) +
  scale_y_continuous(breaks = c(0:3*10, 15)) +
  theme_minimal()

tibble(x = rnorm(150)+1) %>% 
  mutate(y = 5*x+10+rnorm(100, sd = 2)) %>% 
  ggplot(aes(x, y)) +
  geom_point() +
  geom_hline(yintercept = 0, linetype = 2) +
  geom_vline(xintercept = 0, linetype = 2) +
  geom_smooth(method = "lm", se = FALSE) +
  annotate(geom = "point", x = 0, y = 10, size = 4, color = "red") +
  annotate(geom = "label", x = -0.5, y = 12, size = 5, color = "red", label = "intercept") +
  annotate(geom = "label", x = 2, y = 12.5, size = 5, color = "red", label = "slope") +
  scale_x_continuous(breaks = -2:4) +
  scale_y_continuous(breaks = c(0:3*10, 15)) +
  theme_minimal()

## `geom_smooth()` using formula = 'y ~ x'

pBrackets::grid.brackets(815, 420, 815, 545, lwd=2, col="red")

Когда мы пытаемся научиться предсказывать данные одной переменной \(Y\) при помощи другой переменной \(X\), мы получаем формулу:

\[y_i = \hat\beta_0 + \hat\beta_1 \times x_i + \epsilon_i,\] где

• \(x_i\) — \(i\)-ый элемент вектора значений \(X\);
• \(y_i\) — \(i\)-ый элемент вектора значений \(Y\);
• \(\hat\beta_0\) — оценка случайного члена (intercept);
• \(\hat\beta_1\) — оценка углового коэффициента (slope);
• \(\epsilon_i\) — \(i\)-ый остаток, разница между оценкой модели (\(\hat\beta_0 + \hat\beta_1 \times x_i\)) и реальным значением \(y_i\); весь вектор остатков иногда называют случайным шумом (на графике выделены красным).

Причем, иногда мы можем один или другой параметр считать равным нулю.

Определите по графику формулу синей прямой.

Задача регрессии — оценить параметры \(\hat\beta_0\) и \(\hat\beta_1\), если нам известны все значения \(x_i\) и \(y_i\) и мы пытаемся минимизировать значния \(\epsilon_i\). В данном конкретном случае, задачу можно решить аналитически и получить следующие формулы:

\[\hat\beta_1 = \frac{(\sum_{i=1}^n x_i\times y_i)-n\times\bar x \times \bar y}{\sum_{i = 1}^n(x_i-\bar x)^2}\]

\[\hat\beta_0 = \bar y - \hat\beta_1\times\bar x\]

При этом, вне зависимости от статистической школы, у регрессии есть свои ограничения на применение:

• линейность связи между \(x\) и \(y\);
• нормальность распределение остатков \(\epsilon_i\);
• гомоскедастичность — равномерность распределения остатков на всем протяжении \(x\);
• независимость переменных;
• независимость наблюдений друг от друга.

Первая регрессия

количество слов и в рассказах М. Зощенко в зависимости от длины рассказа:

zo <- read_tsv("https://github.com/agricolamz/DS_for_DH/raw/master/data/tidy_zoshenko.csv")
zo %>% 
  filter(word == "и") %>% 
  distinct() %>% 
  ggplot(aes(n_words, n))+
  geom_point()+
  labs(x = "количество слов в рассказе",
       y = "количество и")

Избавимся от выбросов и добавим регрессионную линию при помощи функции geom_smooth():

zo %>% 
  filter(word == "и",
         n_words < 1500) %>% 
  distinct() ->
  zo_filtered
zo_filtered %>%   
  ggplot(aes(n_words, n))+
  geom_point()+
  geom_smooth(method = "lm", se = FALSE)+
  labs(x = "количество слов в рассказе",
       y = "количество и")

Чтобы получить формулу этой линии, нужно запустить функцию, которая оценивает линейную регрессию:

fit0 <- lm(n ~ n_words, data = zo)
fit0

## 
## Call:
## lm(formula = n ~ n_words, data = zo)
## 
## Coefficients:
## (Intercept)      n_words  
##    1.259006     0.006204

fit <- lm(n~n_words, data = zo_filtered)
fit

## 
## Call:
## lm(formula = n ~ n_words, data = zo_filtered)
## 
## Coefficients:
## (Intercept)      n_words  
##    -1.47184      0.04405

\[n_{fit0} = 1.259006 + 0.006204 \times n\_words\]

\[n = -1.47184 + 0.04405 \times n\_words\]

Более подробная информация - в функции summary():

summary(fit0)

## 
## Call:
## lm(formula = n ~ n_words, data = zo)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -20.241  -4.974  -2.995   0.385 133.759 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.259e+00  1.071e-01   11.76   <2e-16 ***
## n_words     6.204e-03  8.304e-05   74.70   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.27 on 46656 degrees of freedom
## Multiple R-squared:  0.1068, Adjusted R-squared:  0.1068 
## F-statistic:  5581 on 1 and 46656 DF,  p-value: < 2.2e-16

summary(fit)

## 
## Call:
## lm(formula = n ~ n_words, data = zo_filtered)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -19.6830  -4.3835   0.8986   4.6486  19.6413 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.471840   2.467149  -0.597    0.553    
## n_words      0.044049   0.003666  12.015   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.945 on 64 degrees of freedom
## Multiple R-squared:  0.6928, Adjusted R-squared:  0.688 
## F-statistic: 144.4 on 1 and 64 DF,  p-value: < 2.2e-16

В разделе Coefficients содержится информация о наших коэффициентах:

• Estimate – полученная оценка коэффициентов;
• Std. Error – стандартная ошибка среднего;
• t value – \(t\)-статистика, полученная при проведении одновыборочного \(t\)-теста, сравнивающего данный коэфициент с 0;
• Pr(>|t|) – полученное \(p\)-значение;
• Multiple R-squared и Adjusted R-squared — одна из оценок модели, показывает связь между переменными. Без поправок совпадает с квадратом коэффициента корреляции Пирсона:

cor(zo_filtered$n_words, zo_filtered$n)^2

## [1] 0.6928376

• F-statistic — \(F\)-статистика полученная при проведении теста, проверяющего, не являются ли хотя бы один из коэффицинтов статистически значимо отличается от нуля. Совпадает с результатами дисперсионного анализа (ANOVA).

Предсказание модели (predict)

Сколько будет “и” в рассказе Зощенко длиной 1000 слов? (на самом деле, такого рассказа нет. Это предсказание)

predict(fit, tibble(n_words = 1000))

##        1 
## 42.57715

Допущения линейной регрессии

Нормальность распределения остатков

Линейная регрессия предполагает нормальность распределения остатков. Когда связь не линейна, то остатки тоже будут распределены не нормально.

Можно смотреть на первый график используя функцию plot() — график остатков.

plot(fit)

Интерпретаций этого графика достаточно много (см. статью про это). * Residuals vs Fitted: тренд в стандартизованных остатках (красная линия) должен идти примерно вдоль пунктирной линии y=0 (нормальность распределения остатков).
* QQplot
* Scale-location plot: показывает гомоскидастичность остатков (постоянство распределения) - красная линия должна идти примерно горизонтально.
* Residuals vs. Leverage: показывает точки наблюдения, наиболее сильно влияющие на коэффициенты модели. Если какая-либо точка вываливается за границы красных пунктирных линий (Cook’s distance), то ее влияние чрезмерно.

Можно смотреть на QQplot:

tibble(res = m1$residuals) %>% 
  ggplot(aes(res))+
  geom_histogram(aes(y = ..density..))+
  stat_function(fun = dnorm, args = list(mean = 0, sd = sd(m1$residuals)), color = "red")
qqnorm(m1$residuals)
qqline(m1$residuals)
tibble(res = m2$residuals) %>% 
  ggplot(aes(res))+
  geom_histogram(aes(y = ..density..))+
  stat_function(fun = dnorm, args = list(mean = 0, sd = sd(m2$residuals)), color = "red")
qqnorm(m2$residuals)
qqline(m2$residuals)
tibble(res = m3$residuals) %>% 
  ggplot(aes(res))+
  geom_histogram(aes(y = ..density..))+
  stat_function(fun = dnorm, args = list(mean = 0, sd = sd(m3$residuals)), color = "red")
qqnorm(m3$residuals)
qqline(m3$residuals)

Гетероскидастичность

Распределение остатков непостоянно (т.е. не гомоскидастично): остатки становятся больше с ростом значений предиктора (или наоборот).

ldt %>% 
  ggplot(aes(Mean_RT, Freq))+
  geom_point()+
  theme_bw()

Тоже решается преобразованием данных.

Мультиколлинеарность

Линейная связь между некоторыми предикторами в модели.

• корреляционная матрица
• VIF (Variance inflation factor), car::vif()
  • VIF = 1 (Not correlated)
  • 1 < VIF < 5 (Moderately correlated)
  • VIF >=5 (Highly correlated) # иногда VIF >=3 или VIF >=10

Независимость наблюдений

Наблюдения должны быть независимы. В ином случае нужно использовать модель со смешанными эффектами.

Линейная модель с категориальными переменными-предикторами

\[y_i = \hat\beta_0 + \hat\beta_1 \times \text{dummy_chekhov} + \epsilon_i,\] \[y_i = \hat\beta_0 + \hat\beta_1 \times 1 + \epsilon_i,\]

dummy_chekhov = 1, если это рассказ Чехова dummy_chekhov = 0, если это рассказ Зощенко (дефолтный)

\[ mean(y_z) = \hat\beta_0 \]

\[ mean(y_{ch}) \]

Давайте рассмотрим простой пример с рассказами Чехова и Зощенко, которые мы рассматривали в прошлом разделе. Мы будем анализировать логарифм доли слов деньги в текстах:

chekhov <- read_tsv("https://github.com/agricolamz/DS_for_DH/raw/master/data/tidy_chekhov.tsv")
zoshenko <- read_tsv("https://github.com/agricolamz/DS_for_DH/raw/master/data/tidy_zoshenko.csv")
chekhov$author <- "Чехов"
zoshenko$author <- "Зощенко"
chekhov_zoshenko <-
  chekhov %>% 
  bind_rows(zoshenko) %>% 
  filter(str_detect(word, "деньг")) %>% 
  group_by(author, titles, n_words) %>% 
  summarise(n = sum(n)) %>% 
  mutate(log_ratio = log(n/n_words)) %>%
  ungroup()

Визуализация выглядит так:

chekhov_zoshenko %>% 
  group_by(author) %>% 
  mutate(mean = mean(log_ratio)) %>% 
  ggplot(aes(author, log_ratio))+
  geom_violin()+
  geom_hline(aes(yintercept = mean), linetype = 3)+
  geom_point(aes(y = mean), color = "red", size = 5)+
  scale_y_continuous(breaks = c(-7, -5, -3, -6.34))

Красной точкой обозначены средние значения, так что мы видим, что между двумя писателями есть разница, но является ли она статистически значимой? В прошлом разделе мы рассмотрели, что в таком случае можно сделать t-test:

t.test(log_ratio~author, 
       data = chekhov_zoshenko, 
       var.equal =TRUE) # поправка Уэлча отключена

## 
##  Two Sample t-test
## 
## data:  log_ratio by author
## t = 5.6871, df = 125, p-value = 8.665e-08
## alternative hypothesis: true difference in means between group Зощенко and group Чехов is not equal to 0
## 95 percent confidence interval:
##  0.8606107 1.7793181
## sample estimates:
## mean in group Зощенко   mean in group Чехов 
##             -5.021262             -6.341226

Разница между группами является статистически значимой (t(125) = 5.6871, p-value = 8.665e-08).

Для того, чтобы запустить регрессию на категориальных данных категориальная переменная автоматически разбивается на группу бинарных dummy-переменных:

tibble(author = c("Чехов", "Зощенко"),
       dummy_chekhov = c(1, 0),
       dummy_zoshenko = c(0, 1))

## # A tibble: 2 × 3
##   author  dummy_chekhov dummy_zoshenko
##   <chr>           <dbl>          <dbl>
## 1 Чехов               1              0
## 2 Зощенко             0              1

Дальше для регрессионного анализа выкидывают одну из переменных, так как иначе модель не сойдется (dummy-переменных всегда n-1, где n — количество категорий в переменной).

tibble(author = c("Чехов", "Зощенко"),
       dummy_chekhov = c(1, 0))

## # A tibble: 2 × 2
##   author  dummy_chekhov
##   <chr>           <dbl>
## 1 Чехов               1
## 2 Зощенко             0

Если переменная dummy_chekhov принимает значение 1, значит речь о рассказе Чехова, а если принимает значение 0, то о рассказе Зощенко. Если вставить нашу переменную в регрессионную формулу получится следующее:

\[y_i = \hat\beta_0 + \hat\beta_1 \times \text{dummy_chekhov} + \epsilon_i,\]

Так как dummy_chekhov принимает либо значение 1, либо значение 0, то получается, что модель предсказывает лишь два значения:

\[y_i = \left\{\begin{array}{ll}\hat\beta_0 + \hat\beta_1 \times 1 + \epsilon_i = \hat\beta_0 + \hat\beta_1 + \epsilon_i\text{, если рассказ Чехова}\\ \hat\beta_0 + \hat\beta_1 \times 0 + \epsilon_i = \hat\beta_0 + \epsilon_i\text{, если рассказ Зощенко} \end{array}\right.\]

Таким образом, получается, что свободный член \(\beta_0\) и угловой коэффициент \(\beta_1\) в регрессии с категориальной переменной получает другую интерпретацию. Одно из значений переменной кодируется при помощи \(\beta_0\), а сумма коэффициентов \(\beta_0+\beta_1\) дают другое значение переменной. Так что \(\beta_1\) – это разница между оценками двух значений переменной.

Запустим регрессию на этих же данных:

fit2 <- lm(log_ratio~author, data = chekhov_zoshenko)
summary(fit2)

## 
## Call:
## lm(formula = log_ratio ~ author, data = chekhov_zoshenko)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.8652 -0.6105 -0.0607  0.6546  3.2398 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -5.0213     0.2120 -23.680  < 2e-16 ***
## authorЧехов  -1.3200     0.2321  -5.687 8.67e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9717 on 125 degrees of freedom
## Multiple R-squared:  0.2056, Adjusted R-squared:  0.1992 
## F-statistic: 32.34 on 1 and 125 DF,  p-value: 8.665e-08

Во-первых, стоит обратить внимание на то, что R сам преобразовал нашу категориальную переменную в dummy-переменную authorЧехов. Во-вторых, можно заметить, что значения t-статистики и p-value совпадают с результатами полученными нами в t-тесте выше. Статистически значимый коэффициент при аргументе authorЧехов следует интерпретировать как разницу средних между логарифмом долей в рассказах Чехова и Зощенко.

Дамми-кодирование с более чем двумя уровнями: \[ y_i = \hat\beta_0 + \hat\beta_1 \times \text{dummy_chekhov} + \hat\beta_2 \times \text{dummy_tolstoy} + \epsilon_i, \] Зощенко Чехова Толстой
dummy_zoschenko 1 0 0
dummy_chekhov 0 1 0
dummy_tolstoy 0 0 1

1 категориальный предиктор (author):
Multiple R-squared: 0.2056, Adjusted R-squared: 0.1992

1 категориальный (author) + 1 непрерывный предиктор (n_word):
Multiple R-squared: 0.4276, Adjusted R-squared: 0.4184

fit3 <- lm(log_ratio ~ author + n_words, data = chekhov_zoshenko)
summary(fit3)

## 
## Call:
## lm(formula = log_ratio ~ author + n_words, data = chekhov_zoshenko)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.5996 -0.6427 -0.0116  0.5475  3.2360 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -4.487e+00  1.964e-01 -22.841  < 2e-16 ***
## authorЧехов -1.088e+00  2.006e-01  -5.422 2.95e-07 ***
## n_words     -6.873e-04  9.909e-05  -6.936 1.97e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8281 on 124 degrees of freedom
## Multiple R-squared:  0.4276, Adjusted R-squared:  0.4184 
## F-statistic: 46.32 on 2 and 124 DF,  p-value: 9.465e-16

Метрика мультиколлинеарности предикторов

car::vif(fit3)

##   author  n_words 
## 1.028633 1.028633

#VIF: 1 / 1-R2

\[R^2 = 1 − Unexplained\_Variation / Total\_Variation\] Примечание: В случае, если один из предикторов - непрерывный, а другой - категориальный (фактор), проверяем количество уровней в факторе. Если уровней всего два, как в данном случае, то vif() интерпретируется обычным образом. Если уровне более двух, то выдаются метрики GVIF and aGSIF (generalized VIF). aGSIF > 1.6 говорит о средней инфляции, aGSIF > 2.2 или aGSIF > 3.2 свидетельствует о большой инфляции.

Вывод модели

# which components are in lm?
names(fit3)

##  [1] "coefficients"  "residuals"     "effects"       "rank"         
##  [5] "fitted.values" "assign"        "qr"            "df.residual"  
##  [9] "contrasts"     "xlevels"       "call"          "terms"        
## [13] "model"

Коэффициенты модели:

fit3$coefficients

##   (Intercept)   authorЧехов       n_words 
## -4.4869205145 -1.0878105849 -0.0006872767

summary(fit3)$coefficients

##                  Estimate   Std. Error    t value     Pr(>|t|)
## (Intercept) -4.4869205145 1.964436e-01 -22.840753 2.951741e-46
## authorЧехов -1.0878105849 2.006117e-01  -5.422467 2.951167e-07
## n_words     -0.0006872767 9.908636e-05  -6.936138 1.971634e-10

Предсказанные значения зависимой переменной:

head(fit3$fitted.values)

##         1         2         3         4         5         6 
## -4.893788 -4.703413 -5.000316 -4.990007 -6.315764 -4.839493

Ошибки модели (остатки, unexplained residuals):

head(fit3$residuals)

##           1           2           3           4           5           6 
## -1.48971833  0.33713439  0.94920037  0.59145111 -1.57069351 -0.01448803

Сумма предсказанных значений и ошибок дает значение зависимой переменной:

chekhov_zoshenko %>%
  select(log_ratio) %>%
  mutate(fit3.predicted = fit3$fitted.values,
         fit3.residuals = fit3$residuals,
         `predicted+residuals` = fit3$fitted.values + fit3$residuals) %>%
  head()

## # A tibble: 6 × 4
##   log_ratio fit3.predicted fit3.residuals `predicted+residuals`
##       <dbl>          <dbl>          <dbl>                 <dbl>
## 1     -6.38          -4.89        -1.49                   -6.38
## 2     -4.37          -4.70         0.337                  -4.37
## 3     -4.05          -5.00         0.949                  -4.05
## 4     -4.40          -4.99         0.591                  -4.40
## 5     -7.89          -6.32        -1.57                   -7.89
## 6     -4.85          -4.84        -0.0145                 -4.85

Дисперсия значений зависимой переменной делится на объясненную моделью и необъясненную. Полная дисперсия (total sum of squares, TSS) может быть подсчитана как сумма квадратов разниц со средним.

tss <- sum((chekhov_zoshenko$log_ratio-mean(chekhov_zoshenko$log_ratio))^2)

Необъясненная дисперсия - сумма квадратов ошибок:

rss <- sum(fit3$residuals^2)

Посчитаем \(R^2\) – долю объясненной дисперсии:

1 - rss/tss

## [1] 0.4276272

#summary(fit3)$r.squared

Модель объясняет 43% дисперсии - много это или мало?

Сравним R2 в моделях fit3 и fit2:

summary(fit2)$r.squared

## [1] 0.2055557

Сравнение моделей

Функция aov() может использоваться для представления выдачи модели в формате выдачи дисперсионного анализа.

aov(fit3)

## Call:
##    aov(formula = fit3)
## 
## Terms:
##                   author  n_words Residuals
## Sum of Squares  30.53837 32.99204  85.03455
## Deg. of Freedom        1        1       124
## 
## Residual standard error: 0.8281078
## Estimated effects may be unbalanced

summary(aov(fit3))

##              Df Sum Sq Mean Sq F value   Pr(>F)    
## author        1  30.54   30.54   44.53 7.45e-10 ***
## n_words       1  32.99   32.99   48.11 1.97e-10 ***
## Residuals   124  85.03    0.69                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Хорошая модель должна обладать свойством парсимонией – не только хорошо обобщать данные, но быть настолько сложной, насколько это необходимо (т. е. не иметь больше предикторов, чем это минимально необходимо). Более сложная модель должна доказать, что она значимо лучше объясняет дисперсию. Если это не так, стоит использовать более простую модель.

anova(fit2,fit3)

## Analysis of Variance Table
## 
## Model 1: log_ratio ~ author
## Model 2: log_ratio ~ author + n_words
##   Res.Df     RSS Df Sum of Sq     F    Pr(>F)    
## 1    125 118.027                                 
## 2    124  85.035  1    32.992 48.11 1.972e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

p-value показывает статистическую значимость добавления новых предикторов. Степень свободы (Df) 1 говорит о том, что две модели различаются на одну переменную.

Множественная регрессия

\[y_i = \hat\beta_0 + \hat\beta_1 \times x_{1i}+ \dots+ \hat\beta_n \times x_{ni} + \epsilon_i,\]

• \(x_{ki}\) — \(i\)-ый элемент векторов значений \(X_1, \dots, X_n\);
• \(y_i\) — \(i\)-ый элемент вектора значений \(Y\);
• \(\hat\beta_0\) — оценка случайного члена (intercept);
• \(\hat\beta_k\) — коэфциент при переменной \(X_{k}\);
• \(\epsilon_i\) — \(i\)-ый остаток, разница между оценкой модели (\(\hat\beta_0 + \hat\beta_1 \times x_i\)) и реальным значением \(y_i\); весь вектор остатков иногда называют случайным шумом.

В такой регресии предикторы могут быть как числовыми, так и категориальными (со всеми вытекающими последствиями, которые мы обсудили в предудщем разделе). Такую регрессию чаще всего сложно визуализировать, так как в одну регрессионную линию вкладываются сразу несколько переменных.

Попробуем предсказать длину лепестка на основе длины чашелистик и вида ириса:

iris %>% 
  ggplot(aes(Sepal.Length, Petal.Length, color = Species))+
  geom_point()

Запустим регрессию:

fit3 <- lm(Petal.Length ~ Sepal.Length+ Species, data = iris)
summary(fit3)

## 
## Call:
## lm(formula = Petal.Length ~ Sepal.Length + Species, data = iris)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.76390 -0.17875  0.00716  0.17461  0.79954 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       -1.70234    0.23013  -7.397 1.01e-11 ***
## Sepal.Length       0.63211    0.04527  13.962  < 2e-16 ***
## Speciesversicolor  2.21014    0.07047  31.362  < 2e-16 ***
## Speciesvirginica   3.09000    0.09123  33.870  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2826 on 146 degrees of freedom
## Multiple R-squared:  0.9749, Adjusted R-squared:  0.9744 
## F-statistic:  1890 on 3 and 146 DF,  p-value: < 2.2e-16

Все предикторы статистически значимы. Давайте посмотрим предсказания модели для всех наблюдений:

iris %>% 
  mutate(prediction = predict(fit3)) %>% 
  ggplot(aes(Sepal.Length, prediction, color = Species))+
  geom_point()

Всегда имеет смысл визуализировать, что нам говорит наша модель. Если использовать пакет ggeffects (или предшествовавший ему пакет effects), это можно сделать не сильно задумываясь, как это делать:

library(ggeffects)
plot(ggpredict(fit3, terms = c("Sepal.Length", "Species")))

Как видно из графиков, наша модель имеет одинаковые угловые коэффициенты (slope) для каждого из видов ириса и разные свободные члены (intercept).

summary(fit3)

## 
## Call:
## lm(formula = Petal.Length ~ Sepal.Length + Species, data = iris)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.76390 -0.17875  0.00716  0.17461  0.79954 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       -1.70234    0.23013  -7.397 1.01e-11 ***
## Sepal.Length       0.63211    0.04527  13.962  < 2e-16 ***
## Speciesversicolor  2.21014    0.07047  31.362  < 2e-16 ***
## Speciesvirginica   3.09000    0.09123  33.870  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2826 on 146 degrees of freedom
## Multiple R-squared:  0.9749, Adjusted R-squared:  0.9744 
## F-statistic:  1890 on 3 and 146 DF,  p-value: < 2.2e-16

\[y_i = \left\{\begin{array}{ll} -1.70234 + 0.63211 \times \text{Sepal.Length} + \epsilon_i\text{, если вид setosa}\\ -1.70234 + 2.2101 + 0.63211 \times \text{Sepal.Length} + \epsilon_i\text{, если вид versicolor} \\ -1.70234 + 3.09 + 0.63211 \times \text{Sepal.Length} + \epsilon_i\text{, если вид virginica} \end{array}\right.\]

Нелинейность взаимосвязи

Давайте восползуемся данными из пакета Rling Натальи Левшиной. В датасете 100 произвольно выбранных слов из проекта English Lexicon Project (Balota et al. 2007), их длина, среднее время реакции и частота в корпусе.

ldt <- read_csv("https://goo.gl/ToxfU6")
ldt

## # A tibble: 100 × 3
##    Length  Freq Mean_RT
##     <dbl> <dbl>   <dbl>
##  1      8   131    819.
##  2     10    82    978.
##  3      7     0    908.
##  4      6   592    766.
##  5     12     2   1125.
##  6     12     9    948.
##  7      3 14013    642.
##  8     11    15    817.
##  9     11    48    931.
## 10      5   290    654.
## # … with 90 more rows

Давайте посмотрим на простой график:

ldt %>% 
  ggplot(aes(Mean_RT, Freq))+
  geom_point()+
  theme_bw()

Регрессия на таких данных будет супер неиформативна:

ldt %>% 
  ggplot(aes(Mean_RT, Freq))+
  geom_point()+
  geom_smooth(method = "lm")+
  theme_bw()

## `geom_smooth()` using formula = 'y ~ x'

m1 <- summary(lm(Mean_RT~Freq, data = ldt))
m1

## 
## Call:
## lm(formula = Mean_RT ~ Freq, data = ldt)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -224.93  -85.42  -30.52   81.90  632.66 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 826.998242  15.229783  54.301  < 2e-16 ***
## Freq         -0.005595   0.001486  -3.765 0.000284 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 143.9 on 98 degrees of freedom
## Multiple R-squared:  0.1264, Adjusted R-squared:  0.1174 
## F-statistic: 14.17 on 1 and 98 DF,  p-value: 0.0002843

Логарифмирование

ldt %>% 
  ggplot(aes(Mean_RT, log(Freq)))+
  geom_point()+
  geom_smooth(method = "lm")+
  theme_bw()

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 5 rows containing non-finite values (`stat_smooth()`).

ldt %>% 
  ggplot(aes(Mean_RT, log(Freq+1)))+
  geom_point()+
  geom_smooth(method = "lm")+
  theme_bw()

## `geom_smooth()` using formula = 'y ~ x'

m2 <- summary(lm(Mean_RT~log(Freq+1), data = ldt))
m2

## 
## Call:
## lm(formula = Mean_RT ~ log(Freq + 1), data = ldt)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -242.36  -76.66  -17.49   48.64  630.49 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    1001.60      29.79  33.627  < 2e-16 ***
## log(Freq + 1)   -34.03       4.76  -7.149 1.58e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 124.8 on 98 degrees of freedom
## Multiple R-squared:  0.3428, Adjusted R-squared:  0.3361 
## F-statistic: 51.11 on 1 and 98 DF,  p-value: 1.576e-10

m1$adj.r.squared

## [1] 0.1174405

m2$adj.r.squared

## [1] 0.336078

Отлогорифмировать можно и другую переменную.

ldt %>% 
  ggplot(aes(log(Mean_RT), log(Freq  + 1)))+
  geom_point()+
  geom_smooth(method = "lm")+
  theme_bw()

## `geom_smooth()` using formula = 'y ~ x'

m3 <- summary(lm(log(Mean_RT)~log(Freq+1), data = ldt))
m1$adj.r.squared

## [1] 0.1174405

m2$adj.r.squared

## [1] 0.336078

m3$adj.r.squared

## [1] 0.3838649

Как интерпретировать полученную регрессию с двумя отлогорифмированными значениями?

В обычной линейной регресии мы узнаем отношения между \(x\) и \(y\): \[y_i = \beta_0+\beta_1\times x_i\]

Как изменится \(y_j\), если мы увеличем \(x_i + 1 = x_j\)? \[y_j = \beta_0+\beta_1\times x_j\]

\[y_j - y_i = \beta_0+\beta_1\times x_j - (\beta_0+\beta_1\times x_i) = \beta_1(x_j - x_i)\]

Т. е. \(y\) увеличится на \(\beta_1\) , если \(x\) увеличится на 1. Что же будет с логарифмированными переменными? Как изменится \(y_j\), если мы увеличем \(x_i + 1 = x_j\)?

\[\log(y_j) - \log(y_i) = \beta_1\times (\log(x_j) - \log(x_i))\]

\[\log\left(\frac{y_j}{y_i}\right) = \beta_1\times \log\left(\frac{x_j}{x_i}\right) = \log\left(\left(\frac{x_j}{x_i}\right) ^ {\beta_1}\right)\]

\[\frac{y_j}{y_i}= \left(\frac{x_j}{x_i}\right) ^ {\beta_1}\]

Т. е. \(y\) увеличится на \(\beta_1\) процентов, если \(x\) увеличится на 1 процент.

Логарифмирование — не единственный вид траснформации:

• трансформация Тьюки

shiny::runGitHub("agricolamz/tukey_transform")

• трансформация Бокса — Кокса
  
• …

Упражнение

В датасет собрана частотность разных лемм на основании корпуса НКРЯ [Ляшевская, Шаров 2009]. Известно, что частотность слова связана с рангом слова (см. закон Ципфа). Постройте переменную ранга и визуализируйте связь ранга и логорифма частотности с разбивкой по частям речи. Какие части речи так и не приобрели после трансформации “приемлемую” линейную форму?

Полезности

• Common statistical tests are linear models link: t-test, ANOVA, correlation explained as lm.

В качестве напоминания

• Дисперсия – мера разброса значений наблюдений относительно среднего.

\[\sigma^2_X = \frac{\sum_{i = 1}^n(x_i - \bar{x})^2}{n - 1},\]

где

• \(x_1, ..., x_n\) — наблюдения;
• \(\bar{x}\) — среднее всех наблюдений;
• \(X\) — вектор всех наблюдений;
• \(n\) — количество наблюдений.

Представим, что у нас есть следующие данные:

set.seed(42)
df <- tibble(x = sort(rnorm(20, mean = 50, sd = 10)), 
             y = seq_along(x))
df %>% 
  ggplot(aes(x)) +
  geom_point(y = 0) +
  ggrepel::geom_text_repel(aes(label = y), y = 0) +
  labs(x = "значение наблюдений x")

Дисперсия - сумма квадратов расстояний от каждой точки до среднего выборки (пунктирная линия) разделенное на n-1.

Распределения могут иметь одинаковое среднее, но разную дисперсию:

x <- rnorm(20, mean = 50, sd = 10)
var(x)

## [1] 107.8731

sd(x) # sqrt(var(x))

## [1] 10.3862

• Ковариация и корреляция - меры ассоциации двух переменных.

\[cov(X, Y) = \frac{\sum_{i = 1}^n(x_i - \bar{x})(y_i-\bar{y})}{n - 1},\] \[\rho_{X,Y} = \frac{cov(X, Y)}{\sigma_X\times\sigma_Y} = \frac{1}{n-1}\times\sum_{i = 1}^n\left(\frac{x_i-\bar{x}}{\sigma_X}\times\frac{y_i-\bar{y}}{\sigma_Y}\right),\]

для коэффициента корреляции Пирсона, где

• \((x_1, y_1), ..., (x_n, y_n)\) — пары наблюдений;
• \(\bar{x}, \bar{y}\) — средние наблюдений;
• \(X, Y\) — векторы всех наблюдений;
• \(n\) — количество наблюдений.

Для таких данных:

Ковариацию можно представить графически следующим образом:

Положительная ковариация - много красных квадратов, отрицательня ковариация - много синих.

Коэффициент корреляции Пирсона можно представить как среднее произведение \(z\)-нормализованных значений двух переменных.

Что можно сказать про корреляцию в следующих данных:

В открытом доступе можно найти игры “Угадай корреляцию” здесь или здесь.